#include "edges.h"
#include "Grid.h"
#include "Parser.h"
#include "hash.h"
#include "Texture2.h"
#include <string>

#define NUMARGS 5

using namespace glm;
using namespace std;

int crash = 0;
int SHADING = 0;
int REFLECT = 6;
int ALIAS = 3;
int INDIRECT = 0;
int ILIMIT = 256;
int LSAMPLE = 25;
int SIM_BOOL = 0;
string SIM = "";

void debug() {
    printf("segfault %d\n", crash++);
}

void genSceneNoise(Scene *s, NoiseGen *n) {
    Geometry *geom;
    
    s->geometry.reset();
    while (s->geometry.hasNext()) {
        geom = s->geometry.getNext();
        geom->setNoise(n);
    }
}

/* for rotating the camera */
mat4 lookTowards(vec3 up, vec3 lookAt, vec3 location) {
    /* centers our coordinates around location */
    vec3 loc = lookAt - location;
    vec3 axis = cross(up, loc);
    mat4 ret = mat4(1.f);
    /* how much to rotate around the Y-axis */
    /* I set it to 0 here because it turned out to be wrong when it wasn't 0 */
    float rotY = (loc.z == 0.0 && loc.x == 0.0) || 1 ? 0.0 :
    atan2(loc.x, loc.z) * 180.0/3.1415926 - 270;
    float rotA = 90 - asin(length(axis)/length(loc)) * 180.0 / 3.1415926;

    if (loc.y > 0) {
        rotA = -rotA;
    }

    ret = rotate(ret, rotA, axis);
    //ret = rotate(ret, rotY, up);
    return ret;
}

inline vec2 jitter(float min, float max) {
    float q;
    vec2 ret;

    q = ((float) rand()) / RAND_MAX;
    ret.x = min - (min - max) * q;
    q = ((float) rand()) / RAND_MAX;
    ret.y = min - (min - max) * q;

    return ret;
}

inline float jitterf(float min, float max) {
    float q;
    float ret;

    q = ((float) rand()) / RAND_MAX;
    ret = min - (min - max) * q;

    return ret;
}

// n is the normal of the geometry we're
// sampling the hemisphere for
vec3 sampleHemisphere(vec3 n, int i) {
    float r, theta;
    float x, y;
    vec3 t, b;

    theta = (2.f * PI * ((float) rand())) / RAND_MAX;
    r = ((float) rand()) / RAND_MAX;

    x = r * cos(theta);
    y = r * sin(theta);
    t = cross(vec3(0, 1, 0), n);
    if (length(t) < 0.1) 
        t = cross(vec3(0, 0, 1), n);

    b = -cross(t, n);

    t = x * t + y * b + sqrt(1.f - r) * n;
    return normalize(t);
}

vec3 refract(Ray *r, Geometry *g) {
    vec3 T, N, I;
    float ndoti, two_ndoti, ndoti2, a, b, b2, D2;

    N = g->getNormal(r);
    I = -r->d;

    ndoti = dot(N, I);
    ndoti2 = ndoti * ndoti;

    if (ndoti >= 0.f) {
        b = g->finish.ior;
    }
    else {
        b = 1.f / g->finish.ior;
        //N = -N;
    }

    b2 = b * b;
    D2 = 1.f - b2 * (1.f - ndoti2);

    if (D2 >= 0.f) {
        if (ndoti >= 0.f) {
            a = b * ndoti - sqrt(D2);
        }
        else {
            a = b * ndoti + sqrt(D2);
        }

        T = a * N - b * I;
    }
    else {
        two_ndoti = ndoti + ndoti;
        T = two_ndoti * N - I;
    }

    return T;
}

class Raytracer {
public: 
    Scene *scene;
    string output;
    Image *image;
    float width, height;
    gNode *planes;
    AABB *headBox;

    Raytracer(int w, int h, string file, Scene *s) {
        Geometry *g;
        AABB *b;
        bNode *boxes;
        int numBoxes = 0;

        boxes = NULL;
        planes = NULL;

        scene = s;
        image = new Image(w, h);
        output = file;
        width = w;
        height = h;

        scene->geometry.reset();
        while (scene->geometry.hasNext()) {
            g = scene->geometry.getNext();
            if (g->type() == T_PLANE) {
                planes = new gNode(planes, (WPlane *) g);
            }
            else if (g->type() == T_CONE) {
            }
            else {
                b = g->generateAABB();
                //g->print();
                //b->print();
                if (b) {
                    boxes = new bNode(boxes, b);
                }
                numBoxes++;
            }
        }

        headBox = generateBoxes(boxes, 0, numBoxes);
    }

    ~Raytracer() {
        gNode *temp;

        /*temp = planes;
        while (temp) {
            delete temp->plane;
        }*/

        delete image;
        //delete scene;
    }

    Geometry *intersectBoxes(AABB *b, Ray *r) {
        static int printNum = 0;
        Geometry *g1, *g2;
        if (b->intersect(r)) {
            if (b->isLeaf) {
                if (b->g->intersect(r)) {
                    return b->g;
                }
                
                return NULL;
            }
            else {
                g1 = intersectBoxes(b->left, r);
                g2 = intersectBoxes(b->right, r);

                if (g2) {
                    return g2;
                } 

                return g1;
            }
        }

/*        b->print();
        printf("Ray:\tpoint: <%.3f, %.3f, %.3f>\tdir <%.3f, %.3f, %.3f>\n",
                r->point.x, r->point.y, r->point.z, r->d.x, r->d.y, r->d.z);
*/
        return NULL;
    }

    Geometry *intersectScene(Ray *r) {
        Geometry *ret = NULL;
        gNode *temp;
        
        /*scene->geometry.reset();
        while(scene->geometry.hasNext()) {
            temp = scene->geometry.getNext();
            if (temp->intersect(r)) {
                ret = temp;
            }
        }*/
        ret = intersectBoxes(headBox, r);

        temp = planes;
        while (temp != NULL) {
            if (temp->plane->intersect(r)) {
                ret = temp->plane;
            }

            temp = temp->next;
        }

        return ret;
    }

    vec3 shadePhong(Geometry *g, WLight *l, Ray *r);
    vec3 shadeGaussian(Geometry *g, WLight *l, Ray *r);
    vec3 shadeBeckmann(Geometry *g, WLight *l, Ray *r);
    vec3 shadeCookTorrance(Geometry *g, WLight *l, Ray *r);
    vec3 shadePhong2(Geometry *g, Geometry *l, Ray *shadow, Ray *r);
    vec3 shadeGaussian2(Geometry *g, Geometry *l, Ray *shadow, Ray *r);
    vec3 shadeBeckmann2(Geometry *g, Geometry *l, Ray *shadow, Ray *r);
    vec3 shadeCookTorrance2(Geometry *g, Geometry *l, Ray *shadow, Ray *r);

    vec3 normalShading(Geometry *g, Ray *r, int ref, int ind, bool o);
    vec3 refractiveShading(Geometry *g, Ray *r, int ref, int ind, bool o);
    vec3 indirectLighting(vec3 point, vec3 normal, vec3 color, int ind);

    vec3 lightRay(Geometry *g, Ray *r, int ref, int ind, bool output) {
        vec3 clr;
        /* beer's law attenuation */
        float atten;
        float dDotd;

        clr = vec3(0.0);

        if (!g) {
            return clr;
        }
        else if (g->finish.lightFinish) {
            return g->pigment.color;
        }
        
        if (g->finish.refraction == 0.0) {
            /* direct color */
            clr = normalShading(g, r, ref, ind, output);

        }
        else {
            clr = refractiveShading(g, r, ref, ind, output);

            /* this ray is entering the material */
            if (0 && dot(-r->d, g->getNormal(r)) < 0) {
                atten = pow(10.f, -r->t * 5000 * 3.345) / 100.f;
                if (atten != 0)
                    clr.x = clr.x * atten;

                atten = pow(10.f, -r->t * 3000 * 3.345) / 100.f;
                if (atten != 0)
                    clr.y = clr.y * atten;

                atten = pow(10.f, -r->t * 20 * 3.345) / 100.f;
                if (atten != 0)
                    clr.z = clr.z * atten;
            }
        }

        if (ind) {
            clr += indirectLighting(
                r->point + r->t * r->d + 0.01f * g->getNormal(r),
                g->getNormal(r), g->pigment.color, ind);
        }

        return clr;
    }


    void trace() {
        Ray r;
        Geometry *g;
        WLight *l;
        vec3 clr, clr_write;
        color_t clr_out;
        int i, j, index;
        float delta, x, y;
        float max;
        float percent = 0.0;
        mat4 lookat;
        vec3 t, u, v, w;
        vec4 tPrime;
        float right = length(scene->camera->right) / 2.0;
        float left = -right;;
        float top = length(scene->camera->up) / 2.0;
        float bot = -top;
        vec2 prev, next;
        vec2 middle;
        clock_t start, cur;
        double diff;

        start = clock();

        lookat = lookTowards(scene->camera->up, scene->camera->lookAt,
                scene->camera->location);
        delta = 1.0 / ALIAS;

        u = normalize(scene->camera->right);
        v = normalize(scene->camera->up);
        w = cross(u, v);

        diff = 0.0;
        
        cout << endl;
        for (i = 0; i < width; i++) {
            for (j = 0; j < height; j++) {
                clr_write = vec3(0);
 
                for (index = 0; index < ALIAS * ALIAS; index++) {
                    x = index % ALIAS;
                    y = index / ALIAS;
                    prev = vec2(x * delta, y * delta);
                    next = vec2(x * delta + delta, y * delta + delta);
                    middle = (prev + next) / 2.0f;
                    if (ALIAS != 1) {
                        middle += jitter(-delta / 2.0, delta / 2.0);
                    }

                    r.point = scene->camera->location;
                    t = r.point + 
                        (left + (right - left) *
                         ((float) i + middle.x) / width) * u +
                        (bot + (top - bot) * 
                         ((float) j + middle.y) / height) * v +
                        -1.0f * w;

                    t = normalize(t - scene->camera->location);
                    tPrime = lookat * vec4(t.x, t.y, t.z, 0.0);
                    t = vec3(tPrime.x, tPrime.y, tPrime.z);

                    r.d = t;
                    r.t = FLT_MAX;

                    g = intersectScene(&r);
                    clr = vec3(0);
                    
                    if (g) {
                        clr = lightRay(g, &r, REFLECT, INDIRECT, true);
                    }

                    /*max = std::max(std::max(clr.x, clr.y), clr.z);
                    if (max > 1.0) {
                        clr /= max;
                    }*/

                    clr_write += clr;
                    cur = clock();
                    if (isatty(1) && 
                            (double) (cur - start) > CLOCKS_PER_SEC / 33.5) {
                        diff += ((double) (cur - start))/CLOCKS_PER_SEC;
                        start = clock();
                        percent = (i * height + j + ((float)index) / 
                                ALIAS * ALIAS) / height / width;
                        printf("Time elapsed: %.2fs\tComplete: %05.2f%%\t\t"
                                "Estimated time remaining: %.2fs\r",
                                diff, percent * 100.0f, diff / percent - diff);
                    }
                }

               clr_write /= ALIAS * ALIAS;

               clr_out.r = clr_write.x;
               clr_out.g = clr_write.y;
               clr_out.b = clr_write.z;

               image->pixel(i, j, clr_out);
            }
        }

        cur = clock();
        diff += ((double) (cur - start)) / CLOCKS_PER_SEC;
        printf("Time elapsed: %.2fs\tComplete: %%100.00\n", diff);
        image->WriteTga((char *)output.c_str());
    }
};

vec3 Raytracer::indirectLighting(vec3 point, vec3 normal, vec3 color, int ind) {
    Ray r;
    vec3 clr;
    int i;

    clr = vec3(0.0);

    for (i = 0; i < ILIMIT / (INDIRECT - ind + 1); i++) {
        r.t = FLT_MAX;
        r.point = point;
        r.d = sampleHemisphere(normal, i);

        clr += lightRay(intersectScene(&r), &r, REFLECT, ind - 1, false);
    }

    clr = color * clr / ((float) ILIMIT);

    return clr;
}

vec3 Raytracer::normalShading(Geometry *g, Ray *r, int ref, int ind, 
        bool output) {
    Ray *rShadow, *rBounce;
    Geometry *gShadow, *gBounce;
    WLight *l;
    Object *lo;
    Geometry *lightObject;
    vec3 clr, clr1, clr2, lclr;
    vec3 normal;
    vec3 loc;
    float dDotn;
    int lightType;
    
    clr = clr1 = clr2 = vec3(0);
    /*if (output) {
        return clr;
    }*/

    rShadow = new Ray();
    scene->lights.reset();
    int count = 0;
    while (scene->lights.hasNext()) {
        lo = scene->lights.getNext();
        lightType = lo->oType;

        rShadow->point = r->point + r->t * r->d;
        rShadow->t = FLT_MAX;

        if (lightType == LIGHT_TYPE) {
            l = (WLight *) lo;

            rShadow->d = 
                normalize(l->location - rShadow->point);
            /* prevent self shadowing */
            rShadow->point += rShadow->d * 0.01f;

            gShadow = intersectScene(rShadow);
            
            if (!gShadow || gShadow->type() == T_PLANE) {
                if (SHADING == 0) {
                    clr = shadePhong(g, l, r);
                }
                else if (SHADING == 1) {
                    clr = shadeGaussian(g, l, r);
                }
                else if (SHADING == 2) {
                    clr = shadeBeckmann(g, l, r);
                }
                else if (SHADING == 3) {
                    clr = shadeCookTorrance(g, l, r);
                }


                clr1 += clr;
            }
        }       
        else {
            lightObject = (Geometry *) lo;

            lclr = vec3(0);
            for (count = 0; count < LSAMPLE; count++) {
                gShadow = NULL;
                loc = lightObject->genRandomPoint(r->point + r->d * r->t,
                                                  count, LSAMPLE);

                rShadow->d =
                    normalize(loc - rShadow->point);
                rShadow->point += rShadow->d * 0.01f;
                gShadow = intersectScene(rShadow);

                if (!gShadow || gShadow->finish.lightFinish) {
                    if (SHADING == 0) {
                        clr = shadePhong2(g, lightObject, rShadow, r);
                    }
                    else if (SHADING == 1) {
                        clr = shadeGaussian2(g, lightObject, rShadow, r);
                    }
                    else if (SHADING == 2) {
                        clr = shadeBeckmann2(g, lightObject, rShadow, r);
                    }
                    else if (SHADING == 3) {
                        clr = shadeCookTorrance2(g, lightObject, rShadow, r);
                    }


                    lclr += clr;
                }
            }

            clr1 += lclr / ((float) LSAMPLE);
        }
    }
    delete rShadow;

    normal = g->getNormal(r);
    dDotn = dot(-r->d, normal);
    if (g->finish.reflection && ref) {
        if (dDotn < 0) {
            normal = -normal;
            dDotn = dot(-r->d, normal);
        }

        rBounce = new Ray();
        rBounce->d = normalize(r->d +
                2.0f * dDotn * normal);
        rBounce->t = FLT_MAX;
        /* move bounce off of the geometry a little */
        rBounce->point = r->point + r->t * r->d +
            normal * 0.01f;

        gBounce = intersectScene(rBounce);
        if (gBounce) {
            clr2 = lightRay(gBounce, rBounce, ref - 1, ind, false);
        }
        delete rBounce;
    }

    dDotn = std::max(std::max(clr1.x, clr1.y), clr1.z);
    if (dDotn > 1) {
        clr1 *= 0.5;
    }


    clr = clr1 + clr2 * g->finish.reflection;
    if (INDIRECT == 0) {
        clr += g->finish.ambient * g->pigment.color;
    }

    return clr;
}

vec3 Raytracer::refractiveShading(Geometry *g, Ray *r, int ref, int ind,
        bool output) {
    Ray *rShadow, *rBounce, *rRefract;
    Geometry *gShadow, *gBounce, *gRefract;
    WLight *l;
    vec3 clr, clr1, clr2, clr3, lclr;
    vec3 normal;
    float deter, fresnel = 1.0;
    float n1, n2;
    float dDotn, R0;


    clr = clr1 = clr2 = clr3 = vec3(0);
/*    if (output) {
        return clr;
    }
*/
    normal = g->getNormal(r);
    dDotn = dot(-r->d, normal);
    /* we're exiting the object */
    if (dDotn < 0) {
        n1 = g->finish.ior;
        n2 = 1.0;
        normal = -normal;
        dDotn = dot(-r->d, normal);
    }
    else {
        n1 = 1.0;
        n2 = g->finish.ior;
    }

    if (ref) {
        rBounce = new Ray();
        rBounce->d = normalize(r->d +
                2.0f * dDotn * normal);
        rBounce->t = FLT_MAX;
        rBounce->point = r->point + r->t * r->d +
                normal * 0.01f;

        gBounce = intersectScene(rBounce);
        if (gBounce) {
            clr1 = lightRay(gBounce, rBounce, ref - 1, ind, false);
        }

        delete rBounce;

        R0 = (n1 - n2) * (n1 - n2) / ((n1 + n2) * (n1 + n2));
        fresnel = R0 + (1.0f - R0) * 
            (1.0f - dDotn) * (1.0f - dDotn) *
            (1.0f - dDotn);
        /*deter = 1.0f - (n1 * n1)/(n2 * n2) * (1.0f - 
                dDotn * dDotn);*/
        /* total internal reflection */
        /*if (deter < 0.f) {
            fresnel = 1.0;
        }
        else {*/
            rRefract = new Ray();
            rRefract->d = refract(r, g);/* normalize(n1/n2 * 
                (r->d + normal * dDotn) * sqrt(1.f + dDotn * dDotn)
                - normal * sqrt(deter));*/

           /* if (rRefract->d.x < 0.001 && rRefract->d.x > -0.001 &&
                rRefract->d.y < 0.001 && rRefract->d.y > -0.001 &&
                rRefract->d.z < 0.001 && rRefract->d.z > -0.001) {
                rRefract->d = r->d;
                clr1 = vec3(1.f, 0.f, 0.f);
            }*/
            rRefract->t = FLT_MAX;
            rRefract->point = r->point + (r->t + 0.01f) * r->d;
        
            gRefract = intersectScene(rRefract);
            if (gRefract) {
                clr2 = lightRay(gRefract, rRefract, ref - 1, ind, false);
            }

            delete rRefract;
       // }
   }

    /*rShadow = new Ray();
    scene->lights.reset();
    int count = 0;
    int lightType;
    Object *lo;
    Geometry *lightObject;
    while (scene->lights.hasNext()) {
        lo = scene->lights.getNext();
        lightType = lo->oType;

        rShadow->point = r->point + r->t * r->d;
        rShadow->t = FLT_MAX;

        if (lightType == LIGHT_TYPE) {
            l = (WLight *) lo;

            rShadow->d = 
                normalize(l->location - rShadow->point);
            /* prevent self shadowing */
      /*      rShadow->point += rShadow->d * 0.01f;

            gShadow = intersectScene(rShadow);
            
            if (!gShadow || gShadow->type() == T_PLANE) {
                if (SHADING == 0) {
                    clr = shadePhong(g, l, r);
                }
                else if (SHADING == 1) {
                    clr = shadeGaussian(g, l, r);
                }
                else if (SHADING == 2) {
                    clr = shadeBeckmann(g, l, r);
                }
                else if (SHADING == 3) {
                    clr = shadeCookTorrance(g, l, r);
                }


                clr3 += clr;
            }
        }       
        else {
            lightObject = (Geometry *) lo;

            lclr = vec3(0);
            for (count = 0; count < LSAMPLE; count++) {
                gShadow = NULL;
                vec3 loc = lightObject->genRandomPoint(r->point + r->d * r->t,
                                                  count, LSAMPLE);

                rShadow->d =
                    normalize(loc - rShadow->point);
                rShadow->point += rShadow->d * 0.01f;
                gShadow = intersectScene(rShadow);

                if (!gShadow || gShadow->finish.lightFinish) {
                    if (SHADING == 0) {
                        clr = shadePhong2(g, lightObject, rShadow, r);
                    }
                    else if (SHADING == 1) {
                        clr = shadeGaussian2(g, lightObject, rShadow, r);
                    }
                    else if (SHADING == 2) {
                        clr = shadeBeckmann2(g, lightObject, rShadow, r);
                    }
                    else if (SHADING == 3) {
                        clr = shadeCookTorrance2(g, lightObject, rShadow, r);
                    }


                    lclr += clr;
                }
            }

            clr3 += lclr / ((float) LSAMPLE);
        }
    }

    delete rShadow;*/

    if (0 && fresnel == 1) {
        clr = vec3(1.f, 0.f, 0.f);
    }
    else {
        clr = clr1 * g->finish.reflection + clr2 * g->finish.refraction + 
            g->pigment.color * (g->finish.diffuse);
    }

    /*if (length(clr) > 1.f)
        clr =  * normalize(clr);*/

    return (clr);
}

vec3 Raytracer::shadePhong(Geometry *g, WLight *l, Ray *r) {
    vec3 normal = normalize(g->getNormal(r));
    vec3 lightDir = normalize(l->location - (r->point + r->t * r->d));
    float NDotL = std::max(0.f, dot(normal, lightDir));
    vec3 R = -lightDir + 2.f * normal * NDotL;
    float VDotR = std::max(0.f, dot(-r->d, R));
    float kspec = 0.0;
    vec3 ret;

    if (g->finish.specular && g->finish.roughness) {
        kspec = pow(VDotR, 1.0f/ g->finish.roughness);
    }

    ret = 
        (g->pigment.color * g->finish.diffuse) * l->color * NDotL +
        (g->pigment.color * g->finish.specular) * l->color * kspec;

//    return vec3(0.0, 0.0, -normal.z);

    return ret;
}

vec3 Raytracer::shadePhong2(Geometry *g, Geometry *l, Ray *shadow, Ray *r) {
    vec3 normal = normalize(g->getNormal(r));
    vec3 lightDir = shadow->d;
    float NDotL = std::max(0.f, dot(normal, lightDir));
    vec3 R = -lightDir + 2.f * normal * NDotL;
    float VDotR = std::max(0.f, dot(-r->d, R));
    float kspec = 0.0;
    vec3 ret;

    if (g->finish.specular && g->finish.roughness) {
        kspec = pow(VDotR, 1.0f/ g->finish.roughness);
    }

    ret = 
        (g->pigment.color * g->finish.diffuse) * l->pigment.color * NDotL +
        (g->pigment.color * g->finish.specular) * l->pigment.color * kspec;

    return ret;
}

vec3 Raytracer::shadeGaussian(Geometry *g, WLight *l, Ray *r) {
    vec3 normal = normalize(g->getNormal(r));
    vec3 lightDir = normalize(l->location - (r->point + r->t * r->d));
    vec3 H = normalize(lightDir - r->d);
    vec3 R = lightDir - 2.f * dot(normal, lightDir) * normal;
    float NDotL = std::max(0.0f, dot(normal, lightDir));
    float NDotH = std::max(0.0f, dot(normal, H));
    float AngleNH = acos(NDotH);
    float VDotR = std::max(0.0f, dot(-r->d, R));
    float kspec = 0.0;
    vec3 ret;
    
    if (g->finish.specular) {
        kspec = exp(-AngleNH * AngleNH /
                    (g->finish.roughness * g->finish.roughness));
    }
 
    ret = 
        (g->pigment.color * g->finish.diffuse) * l->color * NDotL +
        (g->pigment.color * g->finish.specular) * l->color * kspec;

    return ret;
}

vec3 Raytracer::shadeGaussian2(Geometry *g, Geometry *l, Ray *shadow, Ray *r) {
    vec3 normal = normalize(g->getNormal(r));
    vec3 lightDir = shadow->d;
    vec3 H = normalize(lightDir - r->d);
    vec3 R = lightDir - 2.f * dot(normal, lightDir) * normal;
    float NDotL = std::max(0.0f, dot(normal, lightDir));
    float NDotH = std::max(0.0f, dot(normal, H));
    float AngleNH = acos(NDotH);
    float VDotR = std::max(0.0f, dot(-r->d, R));
    float kspec = 0.0;
    vec3 ret;
    
    if (g->finish.specular) {
        kspec = exp(-AngleNH * AngleNH /
                    (g->finish.roughness * g->finish.roughness));
    }
 
    ret = 
        (g->pigment.color * g->finish.diffuse) * l->pigment.color * NDotL +
        (g->pigment.color * g->finish.specular) * l->pigment.color * kspec;

    return ret;
}

vec3 Raytracer::shadeBeckmann(Geometry *g, WLight *l, Ray *r) {
    vec3 normal = normalize(g->getNormal(r));
    vec3 lightDir = normalize(l->location - (r->point + r->t * r->d));
    vec3 H = normalize(lightDir - r->d);
    float NDotL = std::max(0.0f, dot(normal, lightDir));
    float NDotH = std::max(0.0f, dot(normal, H));
    float kspec = 0.0;
    vec3 ret;
    
    if (g->finish.specular && NDotH != 0.0f) {
        kspec = exp(-(1.f-NDotH*NDotH) / 
                (g->finish.roughness * g->finish.roughness *
                 NDotH * NDotH)) / (3.1415926f *
                 g->finish.roughness * g->finish.roughness *
                 NDotH * NDotH * NDotH * NDotH);
    }
 
    ret = 
        (g->pigment.color * g->finish.diffuse) * l->color * NDotL +
        (g->pigment.color * g->finish.specular) * l->color * kspec;

    return ret;
}

vec3 Raytracer::shadeBeckmann2(Geometry *g, Geometry *l, Ray *shadow, Ray *r) {
    vec3 normal = normalize(g->getNormal(r));
    vec3 lightDir = shadow->d;
    vec3 H = normalize(lightDir - r->d);
    float NDotL = std::max(0.0f, dot(normal, lightDir));
    float NDotH = std::max(0.0f, dot(normal, H));
    float kspec = 0.0;
    vec3 ret;
    
    if (g->finish.specular && NDotH != 0.0f) {
        kspec = exp(-(1.f-NDotH*NDotH) / 
                (g->finish.roughness * g->finish.roughness *
                 NDotH * NDotH)) / (3.1415926f *
                 g->finish.roughness * g->finish.roughness *
                 NDotH * NDotH * NDotH * NDotH);
    }
 
    ret = 
        (g->pigment.color * g->finish.diffuse) * l->pigment.color * NDotL +
        (g->pigment.color * g->finish.specular) * l->pigment.color * kspec;

    return ret;
}

vec3 Raytracer::shadeCookTorrance(Geometry *g, WLight *l, Ray *r) {
    vec3 normal = normalize(g->getNormal(r));
    vec3 lightDir = normalize(l->location - (r->point + r->t * r->d));
    vec3 H = normalize(lightDir - r->d);
    float NDotL = std::max(0.0f, dot(normal, lightDir));
    float NDotH = std::max(0.0f, dot(normal, H));
    float AngleNH = acos(NDotH);
    float VDotN = std::max(0.0f, dot(normal, -r->d));
    float VDotH = dot(-r->d, H);
    float beckmann = 0.0;
    float fresnel;
    float gatten = 0.0;
    float kspec = 0.0;
    float R0 = (1.f-g->finish.ior)/(1.f+g->finish.ior);
    vec3 ret;
    
    if (g->finish.specular && NDotH != 0.0f && VDotN != 0.0 && NDotL != 0.0f) {
        if (g->finish.roughness >= 0.1) {
            beckmann = exp(-(1.f-NDotH*NDotH) / 
                (g->finish.roughness * g->finish.roughness *
                 NDotH * NDotH)) / (3.14159265 *
                 g->finish.roughness * g->finish.roughness *
                 NDotH * NDotH * NDotH * NDotH);
        }
        else {
            /* Use gaussian for shiny materials */
            beckmann = exp(-AngleNH * AngleNH /
                    (g->finish.roughness * g->finish.roughness));
        }

        R0 = R0 * R0;
        fresnel = R0 + (1.f - R0) * 
         (NDotH) * (NDotH) * (NDotH) * (NDotH) * (NDotH);

        gatten = std::min(1.f, std::min(2.f*(NDotH)*(VDotN)/(VDotH),
                2.f * (NDotH) * (NDotL) / (VDotH)));

        kspec = beckmann * fresnel * gatten / (3.141592653 * VDotN * NDotL);
    }
 
    ret = 
        g->pigment.color * g->finish.diffuse * l->color * NDotL +
        g->pigment.color * g->finish.specular * l->color * kspec;

    return ret;
}

vec3 Raytracer::shadeCookTorrance2(Geometry *g, Geometry *l, Ray *shadow, 
        Ray *r) {
    vec3 normal = normalize(g->getNormal(r));
    vec3 lightDir = shadow->d;
    vec3 H = normalize(lightDir - r->d);
    float NDotL = std::max(0.0f, dot(normal, lightDir));
    float NDotH = std::max(0.0f, dot(normal, H));
    float AngleNH = acos(NDotH);
    float VDotN = std::max(0.0f, dot(normal, -r->d));
    float VDotH = dot(-r->d, H);
    float beckmann = 0.0;
    float fresnel;
    float gatten = 0.0;
    float kspec = 0.0;
    float R0 = (1.f-g->finish.ior)/(1.f+g->finish.ior);
    vec3 ret;
    
    if (g->finish.specular && NDotH != 0.0f && VDotN != 0.0 && NDotL != 0.0f) {
        if (g->finish.roughness >= 0.1) {
            beckmann = exp(-(1.f-NDotH*NDotH) / 
                (g->finish.roughness * g->finish.roughness *
                 NDotH * NDotH)) / (3.14159265 *
                 g->finish.roughness * g->finish.roughness *
                 NDotH * NDotH * NDotH * NDotH);
        }
        else {
            /* Use gaussian for shiny materials */
            beckmann = exp(-AngleNH * AngleNH /
                    (g->finish.roughness * g->finish.roughness));
        }

        R0 = R0 * R0;
        fresnel = R0 + (1.f - R0) * 
         (NDotH) * (NDotH) * (NDotH) * (NDotH) * (NDotH);

        gatten = std::min(1.f, std::min(2.f*(NDotH)*(VDotN)/(VDotH),
                2.f * (NDotH) * (NDotL) / (VDotH)));

        kspec = beckmann * fresnel * gatten / (3.141592653 * VDotN * NDotL);
    }
 
    ret = 
        g->pigment.color * g->finish.diffuse * l->pigment.color * NDotL +
        g->pigment.color * g->finish.specular * l->pigment.color * kspec;

    return ret;
}

int getIndex(int val, int tx, int ty, int res) {
    switch (val) {
        case 0:
            return tx;
        case 1:
            return tx + (res + 1);
        case 2:
            return tx + 2 * res + 1;
        case 3:
            return tx + res;
        case 4:
            return tx + (res + 1) * (3 * res + 1);
        case 5:
            return tx + (res + 1) * (3 * res + 2);
        case 6:
            return tx + (res + 1) * (3 * res + 2) + res;
        case 7:
            return tx + (res + 1) * (3 * res + 1) + res;
        case 8:
            return ty;
        case 9:
            return ty + 1;
        case 10:
            return ty + res + 1;
        case 11:
            return ty + res + 2;
        default:
            return -1;
    }
}

vec3 getTri(int val, vec3 pos, vec3 delta) {
    vec3 ret;
    switch (val) {
        case 0:
            ret = vec3(pos.x + delta.x / 2.f, pos.y, pos.z);
            break;
        case 1:
            ret = vec3(pos.x + delta.x, pos.y, pos.z + delta.z / 2.f);
            break;
        case 2:
            ret = vec3(pos.x + delta.x / 2.f, pos.y, pos.z + delta.z);
            break;
        case 3:
            ret = vec3(pos.x, pos.y, pos.z + delta.z / 2.f);
            break;
        case 4:
            ret = vec3(pos.x + delta.x / 2.f, pos.y + delta.y, pos.z);
            break;
        case 5:
            ret = vec3(pos.x + delta.x, pos.y + delta.y, pos.z + delta.z / 2.f);
            break;
        case 6:
            ret = vec3(pos.x + delta.x / 2.f, pos.y + delta.y, pos.z + delta.z);
            break;
        case 7:
            ret = vec3(pos.x, pos.y + delta.y, pos.z + delta.z / 2.f);
            break;
        case 8:
            ret = vec3(pos.x, pos.y + delta.y / 2.f, pos.z);
            break;
        case 9:
            ret = vec3(pos.x + delta.x, pos.y + delta.y / 2.f, pos.z);
            break;
        case 10:
            ret = vec3(pos.x + delta.x, pos.y + delta.y / 2.f, pos.z + delta.z);
            break;
        case 11:
            ret = vec3(pos.x, pos.y + delta.y / 2.f, pos.z + delta.z);
            break;
        default:
            ret = pos;
            break;
    }

    return ret;
}

void generateTriangles(Grid *grid, Scene *scene) {
    WTriangle *tri;
    int res;
    GBox *box;
    int x, y, z, i, tri_index;
    vec3 pos;
    int *tri_val;
    int t, tx, ty;
    int count = 0;


    res = grid->res;
    x = y = z = i = 0;
    tx = 0;
    ty = 2 * (res + 1) * (res);

    while (i < res * res * res) {
        box = &grid->boxes[i];
        tri_index = 0;
        tri_val = triTable[box->verts];

        pos = grid->min;

        pos.x += x * grid->delta.x;
        pos.y += y * grid->delta.y;
        pos.z += z * grid->delta.z;


     //   printf("processing: %d out of %d\n", i, res * res * res);

        /*for (int b = 0; b < 16; b++) {
            fprintf(stderr, "%d\n", tri_val[b]);
        }*/
        while (tri_index < 16 && tri_val[tri_index] != -1) {
            int index;
            vec3 val;
            Vertex *vert;

            tri = new WTriangle();
            tri->pigment.color = vec3(0.2588, 1.0, 0.996078);
            tri->finish.diffuse = 0.08f;
            tri->finish.specular = 0.36;
            tri->finish.roughness = 0.020374f;
            tri->finish.reflection = 0.37;
            tri->finish.refraction = 0.58;
            tri->finish.ior = 1.33;
            
            tri->p1 = grid->edge_points[
             grid->getIndex(tri_val[tri_index], tx, ty)].pos;
            tri->p2 = grid->edge_points[
             grid->getIndex(tri_val[tri_index + 1], tx, ty)].pos;
            tri->p3 = grid->edge_points[
             grid->getIndex(tri_val[tri_index + 2], tx, ty)].pos;
            tri->normal = normalize(
             cross(tri->p3 - tri->p1, tri->p2 - tri->p1));
            //tri->surface = true;

            scene->geometry.add(tri);
            tri_index += 3;
        }
        
        i++;
        if (i % (res * res) == 0) {
            x = z = 0;
            y++;
            tx += (res + 1) * (res + 1) + 2 * res + 2;
            ty += (2 * res + 1) * (res + 1) + 2;
        }
        else if (i % res == 0) {
            x = 0;
            z++;
            tx += res + 2;
            ty += 2;
        }
        else {
            x++;
            tx++;
            ty++;
        }
    }

    /*for (i = 0; i < (res + 1) * (res + 1) * res * 3; i++) {
        if (length(grid->verts[i]) != 0.f) {
            grid->verts[i] = normalize(grid->verts[i]);
        }
    }*/

    /*scene->geometry.reset();
    while (scene->geometry.hasNext()) {
        Geometry *o = scene->geometry.getNext();
        WTriangle *t;

        if (!o->name.compare("Triangle")) {
            t = (WTriangle*) o;

            if (t->surface) {
                t->n1 = grid->verts[tri->i1];
                t->n2 = grid->verts[tri->i2];
                t->n3 = grid->verts[tri->i3];
                /*printf("Normal: %f, %f, %f\nn1: %f, %f, %f\n"
                 "n2: %f, %f, %f\nn3: %f, %f, %f\n",
                 t->normal.x, t->normal.y, t->normal.z, 
                 t->n1.x, t->n1.y, t->n1.z, 
                 t->n2.x, t->n2.y, t->n2.z, 
                 t->n3.x, t->n3.y, t->n3.z);*/
 /*           }
        }
    }*/

}

int main(int argc, char *argv[]) {
    string input, temp, line;
    ifstream file;
    int width, height;
    int argOptions = 0;
    Raytracer *tracer;
    Scene *scene;
    NoiseGen *noise;
    ifstream sim_file;
    int frames, particles, i;
    float radius, x, y, z;
    vec3 min, max, cur;
    HashTable *hash;
    HNode head, *htemp;
    Sphere *p;

    /* make sure our random numbers aren't the same every run */
    srand(time(NULL));

    if (argc < NUMARGS) {
        cout << "Usage: raytrace Width Height -I Filename\n";
        exit(1);
    }

    width = atoi(argv[++argOptions]);
    height = atoi(argv[++argOptions]);

    for (argOptions; argOptions < argc; argOptions++) {
        temp = string(argv[argOptions]);
        if (temp.substr(0, 2).compare("-I") == 0) {
        //get the input file name
            if (temp.length() > 2) {
                temp = temp.substr(2);
            }
            else {
                if (argc < argOptions + 1) {
                    cout << "Usage: raytrace Width Height -I Filename\n";
                    exit(1);
                }
                temp = string(argv[++argOptions]);
            }

            input = temp;
        }
        else if (temp.substr(0, 2).compare("-p") == 0) {
            if (temp.length() > 2) {
                temp = temp.substr(2);
            }
            else {
                if (argc < argOptions + 1) {
                    cout << "Usage: raytrace Width Height -I Filename"
                        " [options]\n";
                    exit(1);
                }

                temp = string(argv[++argOptions]);
            }

            SHADING = strtol(temp.c_str(), NULL, 10);
            if (SHADING > 3) {
                cout << "Cannot find shader " << SHADING << endl;
                exit(1);
            }
        }
        else if (temp.substr(0,2).compare("-r") == 0) {
            if (temp.length() > 2) {
                temp = temp.substr(2);
            }
            else {
                if (argc < argOptions + 1) {
                    cout << "Usage: raytrace Width Height -I Filename"
                        " [options]\n";
                    exit(1);
                }

                temp = string(argv[++argOptions]);
            }

            REFLECT = strtol(temp.c_str(), NULL, 10);
        }
        else if (temp.substr(0, 2).compare("-a") == 0) {
            if (temp.length() > 2) {
                temp = temp.substr(2);
            }
            else {
                if (argc < argOptions + 1) {
                    cout << "Usage: raytrace Width Height -I Filename"
                        " [options]\n";
                    exit(1);
                }

                temp = string(argv[++argOptions]);
            }

            ALIAS = strtol(temp.c_str(), NULL, 10);
        }
        else if (temp.substr(0, 2).compare("-n") == 0) {
            if (temp.length() > 2) {
                temp = temp.substr(2);
            }
            else {
                if (argc < argOptions + 1) {
                    cout << "Usage: raytrace Width Height -I Filename"
                        " [options]\n";
                    exit(1);
                }

                temp = string(argv[++argOptions]);
            }

            INDIRECT = strtol(temp.c_str(), NULL, 10);
        }
        else if (temp.substr(0, 2).compare("-c") == 0) {
            if (temp.length() > 2) {
                temp = temp.substr(2);
            }
            else {
                if (argc < argOptions + 1) {
                    cout << "Usage: raytrace Width Height -I Filename"
                        " [options]\n";
                    exit(1);
                }

                temp = string(argv[++argOptions]);
            }

            ILIMIT = strtol(temp.c_str(), NULL, 10);
        }
        else if (temp.substr(0, 2).compare("-w") == 0) {
            if (temp.length() > 2) {
                temp = temp.substr(2);
            }
            else {
                if (argc < argOptions + 1) {
                    cout << "Usage: raytrace Width Height -I Filename"
                        " [options]\n";
                    exit(1);
                }

                temp = string(argv[++argOptions]);
            }

            SIM = temp;
            SIM_BOOL = 1;
        }
    }

    
    if (input.compare("") == 0) {
        cout << "Usage: raytrace Width Height -I Filename\n";
        exit(1);
    }

    if (ALIAS <= 0) {
        ALIAS = 1;
    }

    cout << "Reading the scene\n";
    file.open(input.c_str(), ifstream::in);
    file.unsetf(ios_base::skipws);
    if (file.is_open()) {
        scene = parseScene(file);
    }
    else {
        cout << "Could not open file: " << input << endl;
        exit(1);
    }
    file.close();
    cout << "Read in scene, preparing to trace\n";

    noise = new NoiseGen(256, 256);    
    genSceneNoise(scene, noise);

    if (SIM_BOOL) {
        cout << "Rendering simulation\n";

        sim_file.open(SIM.c_str());
        sim_file >> frames;
        sim_file >> particles;
        sim_file >> radius;

        for (int j = 0; j < frames; j++) {
            head.next = new HNode();
            hash = new HashTable(particles, radius);
            sim_file >> x >> y >> z;
            min = max = vec3(x, y, z);

            p = new Sphere();
            p->center = vec3(x, y, z);
            p->radius = 0.04823f;
            p->id = 0;
            head.next->p = p;
            hash->insert(p);
            p->pigment.color = vec3(0.5, 0.0, 0.4);
            p->pigment.f = 0.3;
        //    scene->geometry.add(p);

            htemp = head.next;

            for (i = 1; i < particles; i++) {
                sim_file >> x >> y >> z;
                if (min.x > x)
                    min.x = x;
                if (min.y > y)
                    min.y = y;
                if (min.z > z)
                    min.z = z;
                if (max.x < x)
                    max.x = x;
                if (max.y < y)
                    max.y = y;
                if (max.z < z)
                    max.z = z;

                htemp->next = new HNode();
                cur = vec3(x, y, z);
                p = new Sphere();
                p->center = cur;
                p->radius = 0.04823f;
                p->id = i;
                htemp->next->p = p;
                hash->insert(p);
                htemp = htemp->next;
                p->pigment.color = vec3(0.5, 0.0, 0.4);
                p->pigment.f = 0.3;
       //         scene->geometry.add(p);
            }

            printf("Making grid\n");
            Grid *grid = new Grid(hash, min, max, 140);
            printf("Generating vertices\n");
            grid->generateVertices();
            printf("Smoothing vertices\n");
            grid->smoothVertices(10);
            //scene->geometry.add(p);
            printf("Generating Triangles\n");
            generateTriangles(grid, scene);
            delete grid;
            stringstream ss;
            ss << j;
            string str = ss.str();
            tracer = new Raytracer(width, height, 
             string(SIM.substr(0, SIM.length() - 3)) + "tga", scene);
            cout << "Tracing " << j << endl;
            tracer->trace();
            cout << "finished\n";
            

            hash->clean();
            delete tracer;
            delete hash;
            file.open(input.c_str(), ifstream::in);
            file.unsetf(ios_base::skipws);
            if (file.is_open()) {
                scene = parseScene(file);
            }
            else {
                cout << "Could not open file: " << input << endl;
                exit(1);
            }
            file.close();
        }
        sim_file.close();
    }
    else {
        //globalScene.debug();
        tracer = new Raytracer(width, height, 
                string(input.substr(0, input.length() - 4)) + ".tga",
                scene);
        cout << "Tracing\n";
        tracer->trace();
        cout << "finished\n";

        delete tracer;
    }

    delete noise;

    return 0;
}
